/* invert a lut
 *
 * Written on: 5/6/01
 * Modified on :
 *
 * 7/7/03 JC
 * 	- generate image rather than doublemask (arrg)
 * 23/3/10
 * 	- gtkdoc
 * 23/5/13
 * 	- fix 1 high input matrices
 * 	- fix file output
 * 4/9/13
 * 	- convert to a class
 */

/*

	This file is part of VIPS.

	VIPS is free software; you can redistribute it and/or modify
	it under the terms of the GNU Lesser General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
	02110-1301  USA

 */

/*

	These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk

 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <glib/gi18n-lib.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <vips/vips.h>

#include "pcreate.h"

/*
#define DEBUG
 */

/* Our state.
 */
typedef struct _VipsInvertlut {
	VipsCreate parent_instance;

	/* Input image.
	 */
	VipsImage *in;

	/* .. and cast to a matrix.
	 */
	VipsImage *mat;

	int size; /* Number of output elements to generate */

	double **data; /* Rows of unpacked matrix */
	double *buf;   /* Output buffer */
} VipsInvertlut;

typedef VipsCreateClass VipsInvertlutClass;

G_DEFINE_TYPE(VipsInvertlut, vips_invertlut, VIPS_TYPE_CREATE);

static void
vips_invertlut_dispose(GObject *gobject)
{
	VipsInvertlut *lut = (VipsInvertlut *) gobject;

	VIPS_FREE(lut->data);
	VIPS_FREE(lut->buf);
	VIPS_UNREF(lut->mat);

	G_OBJECT_CLASS(vips_invertlut_parent_class)->dispose(gobject);
}

/* Use this to sort our input rows by the first column.
 */
static int
compare(const void *a, const void *b)
{
	double **r1 = (double **) a;
	double **r2 = (double **) b;

	double diff = r1[0][0] - r2[0][0];

	if (diff > 0)
		return 1;
	else if (diff == 0)
		return 0;
	else
		return -1;
}

static int
vips_invertlut_build_init(VipsInvertlut *lut)
{
	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(lut);

	int x, y;

	if (!lut->mat ||
		lut->mat->Xsize < 2 ||
		lut->mat->Ysize < 1) {
		vips_error(class->nickname, "%s", _("bad input matrix"));
		return -1;
	}
	if (lut->size < 1 ||
		lut->size > 65536) {
		vips_error(class->nickname, "%s", _("bad size"));
		return -1;
	}

	if (!(lut->buf =
				VIPS_ARRAY(NULL, lut->size * (lut->mat->Xsize - 1), double)))
		return -1;

	if (!(lut->data = VIPS_ARRAY(NULL, lut->mat->Ysize, double *)))
		return -1;
	for (y = 0; y < lut->mat->Ysize; y++)
		lut->data[y] = VIPS_MATRIX(lut->mat, 0, y);

	/* Sanity check for data range.
	 */
	for (y = 0; y < lut->mat->Ysize; y++)
		for (x = 0; x < lut->mat->Xsize; x++)
			if (lut->data[y][x] > 1.0 ||
				lut->data[y][x] < 0.0) {
				vips_error(class->nickname,
					_("element (%d, %d) is %g, outside range [0,1]"),
					x, y, lut->data[y][x]);
				return -1;
			}

	/* Sort by 1st column in input.
	 */
	qsort(lut->data, lut->mat->Ysize, sizeof(double *), compare);

#ifdef DEBUG
	printf("Input table, sorted by 1st column\n");
	for (y = 0; y < lut->mat->Ysize; y++) {
		printf("%.4d ", y);

		for (x = 0; x < lut->mat->Xsize; x++)
			printf("%.9f ", lut->data[y][x]);

		printf("\n");
	}
#endif /*DEBUG*/

	return 0;
}

static int
vips_invertlut_build_create(VipsInvertlut *lut)
{
	int bands = lut->mat->Xsize - 1;
	int height = lut->mat->Ysize;

	int b;

	/* Do each output channel separately.
	 */
	for (b = 0; b < bands; b++) {
		/* The first and last lut positions we know real values for.
		 */
		int first = lut->data[0][b + 1] * (lut->size - 1);
		int last = lut->data[height - 1][b + 1] * (lut->size - 1);

		int k;

		/* Extrapolate bottom and top segments to (0,0) and (1,1).
		 */
		for (k = 0; k < first; k++) {
			/* Have this inside the loop to avoid /0 errors if
			 * first == 0.
			 */
			double fac = lut->data[0][0] / first;

			lut->buf[b + k * bands] = k * fac;
		}

		for (k = last; k < lut->size; k++) {
			/* Inside the loop to avoid /0 errors for last ==
			 * (size - 1).
			 */
			double fac = (1 - lut->data[height - 1][0]) /
				((lut->size - 1) - last);

			lut->buf[b + k * bands] =
				lut->data[height - 1][0] + (k - last) * fac;
		}

		/* Interpolate the data sections.
		 */
		for (k = first; k <= last; k++) {
			/* Where we're at in the [0,1] range.
			 */
			double ki = (double) k / (lut->size - 1);

			double irange, orange;
			int j;

			/* Search for the lowest real value < ki. There may
			 * not be one: if not, just use 0. Tiny error.
			 */
			for (j = height - 1; j >= 0; j--)
				if (lut->data[j][b + 1] < ki)
					break;
			if (j == -1)
				j = 0;

			/* Interpolate k as being between row data[j] and row
			 * data[j + 1].
			 */
			irange = lut->data[j + 1][b + 1] - lut->data[j][b + 1];
			orange = lut->data[j + 1][0] - lut->data[j][0];

			lut->buf[b + k * bands] = lut->data[j][0] +
				orange * ((ki - lut->data[j][b + 1]) / irange);
		}
	}

	return 0;
}

static int
vips_invertlut_build(VipsObject *object)
{
	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS(object);
	VipsCreate *create = VIPS_CREATE(object);
	VipsInvertlut *lut = (VipsInvertlut *) object;

	if (VIPS_OBJECT_CLASS(vips_invertlut_parent_class)->build(object))
		return -1;

	if (vips_check_matrix(class->nickname, lut->in, &lut->mat))
		return -1;

	if (vips_invertlut_build_init(lut) ||
		vips_invertlut_build_create(lut))
		return -1;

	vips_image_init_fields(create->out,
		lut->size, 1, lut->mat->Xsize - 1,
		VIPS_FORMAT_DOUBLE, VIPS_CODING_NONE,
		VIPS_INTERPRETATION_HISTOGRAM, 1.0, 1.0);
	if (vips_image_write_line(create->out, 0, (VipsPel *) lut->buf))
		return -1;

	return 0;
}

static void
vips_invertlut_class_init(VipsInvertlutClass *class)
{
	GObjectClass *gobject_class = G_OBJECT_CLASS(class);
	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS(class);

	gobject_class->dispose = vips_invertlut_dispose;
	gobject_class->set_property = vips_object_set_property;
	gobject_class->get_property = vips_object_get_property;

	vobject_class->nickname = "invertlut";
	vobject_class->description = _("build an inverted look-up table");
	vobject_class->build = vips_invertlut_build;

	VIPS_ARG_IMAGE(class, "in", 0,
		_("Input"),
		_("Matrix of XY coordinates"),
		VIPS_ARGUMENT_REQUIRED_INPUT,
		G_STRUCT_OFFSET(VipsInvertlut, in));

	VIPS_ARG_INT(class, "size", 5,
		_("Size"),
		_("LUT size to generate"),
		VIPS_ARGUMENT_OPTIONAL_INPUT,
		G_STRUCT_OFFSET(VipsInvertlut, size),
		1, 1000000, 256);
}

static void
vips_invertlut_init(VipsInvertlut *lut)
{
	lut->size = 256;
}

/**
 * vips_invertlut: (method)
 * @in: input mask
 * @out: (out): output LUT
 * @...: `NULL`-terminated list of optional named arguments
 *
 * Given a mask of target values and real values, generate a LUT which
 * will map reals to targets.
 *
 * Handy for linearising images from measurements of a colour chart. All
 * values in [0,1]. Piecewise linear interpolation, extrapolate head and tail
 * to 0 and 1.
 *
 * Eg. input like this:
 *
 * ```
 * 4 3
 * 0.1 0.2 0.3 0.1
 * 0.2 0.4 0.4 0.2
 * 0.7 0.5 0.6 0.3
 * ```
 *
 * Means a patch with 10% reflectance produces an image with 20% in
 * channel 1, 30% in channel 2, and 10% in channel 3, and so on.
 *
 * Inputs don't need to be sorted (we do that). Generate any precision
 * LUT, default to 256 elements.
 *
 * It won't work too well for non-monotonic camera responses
 * (we should fix this). Interpolation is simple piecewise linear; we ought to
 * do something better really.
 *
 * ::: tip "Optional arguments"
 *     * @size: `gint`, generate this much
 *
 * ::: seealso
 *     [method@Image.buildlut].
 *
 * Returns: 0 on success, -1 on error
 */
int
vips_invertlut(VipsImage *in, VipsImage **out, ...)
{
	va_list ap;
	int result;

	va_start(ap, out);
	result = vips_call_split("invertlut", ap, in, out);
	va_end(ap);

	return result;
}
